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A  general  purpose  3D  finite  element  method  model  has  been  developed  for  the  estimation  of  the  compres¬ 
sion  pressure  distribution  in  fuel  cell  stacks.  The  model  can  be  used  for  the  optimisation  of  any  type  of  fuel 
cell  structure  at  any  temperature.  The  model  was  validated  with  pressure  sensitive  film  measurements 
using  PEFC  stack  components  that  had  low  rigidity  and  were  highly  deformable. 
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1.  Introduction 

Uniform  compression  pressure  distribution  over  the  active  area 
and  gaskets  of  the  fuel  cell  stack  is  essential  to  its  performance  and 
life-time  [1]  as  it  has  a  major  impact  on  the  properties  of  the  gas 
diffusion  layer  (GDL)  and  contact  resistances  [2-10],  affecting  both 
current  and  temperature  distributions.  The  compression  pressure 
is  uneven  both  under  the  ribs  in  the  flow  fields  [1,11-13],  as  well  as 
over  the  whole  membrane  electrode  assembly  (MEA)  and  gaskets 
[14-16]. 

For  the  low  temperature  fuel  cells  the  pressure  distribution 
can  be  measured  experimentally  near  operating  temperature  using 
pressure  sensitive  films.  For  high  temperature  fuel  cells  this  is 
not  possible,  due  to  a  limited  operating  temperature  range  of  the 
pressure  sensitive  films.  Therefore,  the  compression  pressure  distri¬ 
bution  within  a  high  temperature  fuel  cell  stack  must  be  modelled, 
which  creates  a  need  for  an  experimentally  validated  modelling 
tool.  With  a  high  quality  modelling  tool,  optimal  end  plate  struc¬ 
tures  and  clamping  schemes  of  both  high  and  low  temperature  fuel 
cells  can  be  readily  designed. 

In  the  past  a  few  attempts  have  been  made  to  simulate  the  com¬ 
pression  pressure  within  a  polymer  electrolyte  fuel  cell  (PEFC)  using 
finite  element  method  (FEM)  computer  models  [14,15,21  ].  Lee  et  al. 


Abbreviations:  EQL,  equalization  layer;  FEM,  finite  element  method;  GDL,  gas 
diffusion  layer;  MEA,  membrane  electrode  assembly;  PEFC,  polymer  electrolyte  fuel 
cell. 
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[14]  modelled  the  pressure  distribution  inside  a  small  unit  cell  and 
validated  their  model  with  experiments.  However,  there  was  up 
to  a  60%  difference  between  model  predictions  and  measurement 
results,  but  quantitatively,  the  trends  in  results  were  similar.  Liu 
et  al.  [15]  presented  a  method  to  optimise  the  assembly  pressure 
and  clamping  bolt  locations  to  achieve  more  uniform  pressure  dis¬ 
tribution  in  a  unit  cell.  Their  work  did  not  contain  experimental 
validation.  Karvonen  et  al.  used  a  FEM  model  to  investigate  pres¬ 
sure  distribution  in  a  multicell  stack,  and  compared  various  end 
plate  structures  and  materials  and  clamping  force  schemes  [17]. 
Their  validation  process  did  not  include  numerical  comparison,  but 
qualitatively  the  results  were  in  agreement.  The  developed  models 
have  shown  the  importance  of  correct  end  plate  design  and  clamp¬ 
ing  pressure  application,  but  none  of  them  can  reliably  predict  the 
compression  pressure  distribution  inside  a  cell  or  a  stack. 

The  purpose  of  this  work  was  to  develop  and  validate  a  widely 
applicable  model  that  can  accurately  predict  pressure  distribu¬ 
tions  within  any  type  of  fuel  cell  at  any  temperature.  A  3D  FEM 
computer  model  was  developed  and  validated  through  experi¬ 
ments  with  pressure  sensitive  films  using  in-house  PEFC  stack  with 
highly  deformable  components.  Results  from  validation  experi¬ 
ments  proved  the  applicability  of  the  model. 

2.  Description  of  the  model 

2.1.  Theory 

The  equation  which  is  solved  in  the  model  is: 

V  •  (DVu)  =  0  (1) 
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Nomenclature 

D  elasticity  matrix 

E  Young’s  modulus  (Pa) 

u  deformation  vector 

u ,  v}  w  deformation  vector  components  (m) 
a  stress  (Nm-2) 

v  Poisson’s  ratio 


where  D  is  the  elasticity  matrix  and  u  is  the  deformation  vector 
consisting  of  deformations  u ,  v,  and  w  in  x-,  y-  and  z-directions, 
respectively: 

u  =  ui  +  vj  +  wk  (2) 


From  the  solution,  stress  and  strain  components  can  be  calculated. 
In  this  work,  we  are  mostly  interested  in  the  perpendicular  stress 
on  stack  components.  For  example,  stress  a  in  the  z-direction  for 
isotropic  materials  can  be  written  as: 


o-z 


E 

(1  +v)(2-v) 


dv 

+  9y  +  1 


(3) 


where  E  is  Young’s  modulus  and  v  Poisson’s  ratio.  In  the  case  of 
anisotropic  materials,  Young’s  modulus  and  Poisson’s  ratio  are  not 
constants  and  thus  the  elasticity  matrix  differs  from  that  of  the 


isotropic  case.  More  information  on  the  subject  can  be  found  in 
reference  [17]  or  structural  mechanics  text  books,  for  example  ref. 
[18]. 

2.2.  Model  description 

The  fuel  cell  model  geometry  used  in  this  model  is  based  on  an 
existing  fuel  cell  stack  design  on  which  experiments  could  readily 
be  carried  out.  The  stack  structure  was  simplified  for  computa¬ 
tional  and  experimental  validation  purposes  by  omitting  bipolar 
plates,  porous  transport  layers  and  gaskets  and  including  only 
end  plates,  flexible  equalization  layers  (EQLs),  current  collectors 
and  bipolar  plate  substitute.  The  resulting  structure  is  artificial 
but  still  corresponds  closely  to  the  situation  in  the  first  and  last 
cells  in  the  real  stacks,  where  the  compression  pressure  distri¬ 
bution  is  usually  the  most  uneven.  Further  into  the  stack,  the 
flexibility  of  cell  components  makes  the  pressure  distribution  more 
even. 

The  model  is  comprised  of  five  separate  domains,  which  were 
the  40  mm  thick  end  plate,  5  mm  thick  equalization  layer,  10  mm 
current  collector,  1  mm  bipolar  plate  substitute  and  gas  line  con¬ 
nector  substitutes.  The  EQL  was  optional,  and  was  not  used  in  all 
the  modelled  cases.  The  model  geometry  is  shown  in  Fig.  1.  The 
geometry  of  each  component  is  identical  to  the  real  geometry  of  the 
actual  stack  components  except  the  bipolar  plate  substitute,  which 
is  identical  in  x-y  dimension,  but  thinner  and  does  not  contain  flow 
channels. 


Fig.  1.  Left:  Stack  overview.  Right:  Stack  quarter  section  view  along  symmetry  planes. 


Table  1 

Material  parameters  used  in  the  model. 


Component/dimensions 

Material 

Young’s  modulus  (Pa) 

Poisson’s  ratio 

Sources 

End  plate  (130  mm  x  96  mm  x  40  mm) 

Ensigner  TECATRON 

GF40  PPS  +  40%  glass 
fibre 

13  x  109 

0.36 

[19] 

Current  collector  (127  mm  x  93  mm  x  10  mm) 

GrafTech  GrafCell 
FFP-300 

1.4  x  109 

0.25 

[20] 

EQL  (115  mm  x  82  mm  x  5  mm) 

Butyl  rubber 

7.2  x  106 

0.49 

Young’s  modulus  measured3,  Poisson’s 
ratio  estimated 

Bipolar  plate  substitute 
(115  mm  x  81  mm  x  1  mm) 

PTFE 

0.4  x  109 

0.25 

COMSOL  material  library 

Gas  line  fitting  substitutes  cylindrical, 
r\  =  10  mm,  r2  =  7  mm,  h  =  40  mm 

PTFE 

0.4  x  109 

0.25 

COMSOL  material  library 

a  Measured  for  0-12  bar  compression  using  Testometric  M500-50  CT  universal  material  tester.  The  response  in  that  range  is  linear.  The  same  value  has  been  used  for  all 
compression  pressures. 
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Fig.  2.  Left:  Model  geometry  and  symmetry  planes  (1)  End  plate  (2)  Equalization  layer  (3)  Current  collector  (4)  Bipolar  plate  substitute  (5)  Gas  line  connector  substitute. 
Right:  Areas  used  for  clamping  force  boundary  conditions. 


Due  to  symmetry,  the  model  implements  only  one  eighth  of  the 
stack.  Symmetry  conditions  were  applied  onto  all  internal  section 
boundaries.  A  symmetry  boundary  condition  was  also  applied  onto 
the  bottom  of  the  bipolar  plate  substitute,  corresponding  to  the 
experimental  scenario  where  the  stack  is  assembled  with  a  thin 
layer  of  PTFE  in  place  of  bipolar  plates.  In  the  model,  clamping  forces 
were  applied  on  a  washer-shaped  area  around  each  bolt  hole.  The 
inner  and  outer  diameters  were  9  and  16  mm,  respectively,  giving 
an  area  of  1.37  x  10-4  m2.  All  other  boundaries  were  free. 

Young’s  modulus  and  Poisson’s  ratio  were  needed  for  each  com¬ 
ponent  to  model  their  deformation  under  stress.  Manufacturer 
supplied  and  literature  data  was  used  whenever  available.  Table  1 
contains  a  list  of  modelled  components,  their  dimensions,  materi¬ 
als  and  material  parameters  at  room  temperature,  and  data  sources. 
All  components  were  assumed  isothermal,  elastic  and  isotropic.  For 
the  end  plate,  EQL  and  bipolar  plate  substitute  isotropy  is  a  reason¬ 
ably  good  assumption,  but  the  properties  of  the  current  collector 
material  are  clearly  different  in  through-  and  in-plane  directions. 
However,  the  error  caused  by  the  assumption  should  be  insignifi¬ 
cant.  The  deformation  of  the  current  collector  is  minor  compared 
to  the  EQL  and  thus,  current  collector’s  properties  in  the  in-plane 
direction  have  negligible  contribution  to  the  resulting  pressure  dis¬ 
tribution. 

The  model  was  implemented  using  COMSOL  Multiphysics  3.4, 
a  commercial  FEM  solver,  running  on  a  computer  with  a  dual  core 
AMD  Opteron  265  processor,  16  GB  of  RAM,  and  GNU/Linux  Open¬ 
SUSE  10.3  64-bit  operating  system. 


3.  Experimental 

3.1.  Description  of  the  pressure  distribution  measurement 
hardware 

For  the  validation  measurements,  a  new  metal-free  end  plate 
structure  was  developed  for  a  previously  designed  in-house  stack. 
The  structure  included  relatively  soft  end  plates  made  of  fibre 
glass  reinforce  polyphenylene  sulphide  (PPS,  Ensigner  Tecatron  GF 
40),  expanded  natural  graphite  current  collectors  (GrafTech  Graf- 
Cell  FFP-300)  and  butyl  rubber  (polyisobutadiene)  EQLs.  The  EQLs 
between  end  plates  and  current  collectors  deform  and  thereby 
transfer  pressure  from  the  high  pressure  regions  to  the  low  pressure 
regions.  The  developed  new  end  plate  structure  is  metal  free.  The 
motivation  for  the  development  of  the  metal-free  structure  comes 
from  the  fact  that  cations,  especially  Fe2+  and  Cu2+,  released  from 
any  other  metal  part  of  the  system  can  act  as  catalysts  for  radi¬ 
cal  formation  from  H202  and  also  decrease  the  conductivity  of  the 


membrane  [21-24].  Therefore,  to  maximise  polymer  fuel  cell  life¬ 
time  it  is  desirable  to  eliminate  all  metallic  components  in  parts  of 
the  fuel  cell  system  that  are  in  contact  with  humidified  reactant  gas 
flows. 

The  experimental  set-up  is  shown  in  Fig.  2.  Between  the  current 
collectors  a  bipolar  plate  substitute  (two  0.5  mm  thick  PTFE  sheets) 
was  added  to  reduce  the  disturbing  effect  of  the  current  collec¬ 
tor  roughness  on  the  pressure  distribution  measurements  and  to 
avoid  discolouration  of  the  pressure  sensitive  films  by  the  current 
collectors.  Pressure  sensitive  films  (Pressurex®  Super  Low  Pressure 
LLW,  FUJIFILM  NDT  Systems,  Inc.)  were  placed  within  the  bipolar 
plate  substitute  (see  Fig.  2).  The  compression  pressure  was  exerted 
by  a  stack  of  disc  springs  (Mubea  Tellerfedern  GmbH)  around  each 
clamping  bolt. 

A  set  of  calibration  measurements,  in  the  range  of  3-28  bar,  were 
performed  to  obtain  a  relation  between  the  colour  intensity  in  the 
pressure  sensitive  film  and  applied  pressure.  The  resulting  calibra¬ 
tion  curve  was  used  in  data  processing  to  convert  the  colour  pattern 
on  the  pressure  sensitive  film  to  numeric  pressure  distribution  data. 
The  applicable  pressure  range  was  found  to  be  4-25  bar. 

All  experimental  work  reported  here  was  carried  out  under  con¬ 
trolled  ambient  conditions.  Air  temperature  and  relative  humidity 
were  maintained  at  23  °C  and  50%,  respectively.  Impact  time  for  all 
pressure  sensitive  film  measurements  was  2  h.  All  measurements 
were  performed  twice. 

3.2.  Test  procedure  for  validation  of  the  computer  model 

Five  test  cases  (see  Table  2)  were  compiled  for  model  valida¬ 
tion.  Pressure  distribution  between  the  bipolar  plate  substitutes 
was  both  modelled  and  experimentally  measured  for  all  cases.  The 
effect  of  equalization  layer  was  investigated  with  three  test  cases: 
no  equalization  layer,  basic  equalization  layer  and  modified  equal¬ 
ization  layer.  The  basic  EQL  is  solid,  except  for  the  holes  required  for 
clamping  bolts  and  gas  channels.  The  modified,  but  non-optimised 
EQL  contains  additional  holes,  which  modify  the  overall  rigidity  in 
the  vicinity  of  the  holes  (see  Fig.  3).  Furthermore,  Fig.  3  shows  the 


Table  2 

Test  cases  for  model  validation. 


EQL 

Clamping  force  arrangement 

3 1<N  at  each  bolt 

High  forces  on 
short  edge 

High  forces  on 
long  edge 

No  EQL 

Case  1 

Basic 

Case  2 

Modified 

Case  3 

Case  4 

Case  5 

272 


M.  Mikkola  et  al.  /  Journal  of  Power  Sources  193  (2009)  269-275 


Fig.  3.  Left:  clamping  force  cases.  Top:  Case  4  and  bottom:  Case  5.  Right:  equalization  layer.  The  basic  EQL  does  not  have  the  perforations  marked  with  darker  colour.  Below: 
the  locations  of  grid  squares  are  given  using  the  coordinate  system  shown  in  this  figure. 


Fig.  4.  (a-e)  Cases  1-5,  respectively.  Modelled  pressure  distributions  on  the  bottom  of  the  bipolar  plate  substitute  over  a  90  mm  x  60  mm  rectangle  on  the  active  area.  The 
centre  of  the  stack  is  at  the  lower  left  corner  of  each  picture. 
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location  of  6  x  4  grid  of  15  mm  by  15  mm  squares,  into  which  the 
active  area  was  divided.  This  grid  was  used  in  the  model  validation 
process. 

Three  test  cases  were  also  used  for  clamping  force  distribution. 
Case  3  is  a  base  case  in  which  3  kN  force  was  applied  to  all  clamping 
bolts.  In  Cases  4  and  5,  the  majority  of  the  force  was  applied  around 
the  bolts  on  the  short  and  long  edges  of  the  stack,  respectively. 
Clamping  forces  and  their  locations  can  be  seen  in  Fig.  3.  In  all 
three  cases  the  total  force  applied  onto  one  quarter  of  the  end  plate 
is  12  kN.  This  corresponds  to  an  average  pressure  of  12.5  bar  on  the 
EQL  and  bipolar  plate  substitute,  which  both  have  an  area  of  ca. 
96  cm2. 

The  results  from  the  model  and  the  measured  data  were  ana¬ 
lyzed  in  the  following  way.  Resulting  pressure  maps  were  scanned 
at  300  DPI  to  digital  format  and  processed  using  Matlab®.  Since 
the  system  is  symmetric  the  measured  data  was  averaged  over  one 
quarter  of  the  active  area  in  order  to  get  a  larger  statistical  selection. 
To  reduce  noise,  the  data  from  the  scans  was  further  averaged  over 
squares  so  that  the  final  images  had  ca.  104  data  points 

The  colouration  of  the  pressure  maps  were  converted  into  sur¬ 
face  pressures  using  the  obtained  calibration  curve.  Since  data 


measured  above  the  saturation  limit  of  the  pressure  sensitive  film 
cannot  be  considered  reliable  all  data  indicating  pressures  higher 
than  the  saturation  limit  were  set  to  25  bar.  To  analyze  the  results, 
the  compression  pressure  over  the  active  area  was  integrated  over 
the  grid  squares  shown  in  Fig.  3,  and  resulting  forces  were  compared 
between  the  measurement  and  modelling  results. 

4.  Results  and  discussion 

4.1.  Modelling  results 

Modelled  clamping  pressure  distributions  over  the  active  area 
are  presented  in  Fig.  4  for  different  cases  given  in  Table  2.  Although 
the  pressure  distribution  was  modelled  and  measured  over  the 
whole  area  of  the  bipolar  plate  substitute,  only  the  distribution  on 
the  active  area  is  shown  in  the  figures.  In  Case  1,  there  is  no  EQL 
present  the  system  to  transfer  pressure  from  high  pressure  regions 
to  low  pressure  regions  and  the  compression  pressure  is  applied  to 
the  outer  edges  of  the  active  area,  especially  the  outer  corner.  In 
Case  2,  a  solid  EQL  is  applied  which  compensates  for  the  deforma¬ 
tion  of  the  end  plate  and  transfers  some  pressure  from  the  corners 


Fig.  5.  (a-e)  Cases  1-5,  respectively.  Measured  pressure  distribution  on  the  bottom  of  the  BPP  substitute  over  a  90  mm  x  60  mm  rectangle  on  the  active  area.  The  centre  of 
the  stack  is  at  the  lower  left  corner  of  each  picture.  Top:  varying  EQL,  bottom:  varying  clamping  force  scheme.  Note  the  scale  in  each  subpicture. 
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towards  the  middle  of  the  active  area.  In  Case  3,  an  EQL,  which  is 
modified  for  increased  deformation,  distributes  the  pressure  even 
more  evenly  over  the  active  area. 

The  compression  pressure  distribution  caused  by  applying  dif¬ 
ferent  forces  to  the  bolts  in  the  stack  was  studied  with  modified  EQL. 
Comparing  test  Cases  3-5,  it  can  be  seen  that  the  clamping  force 
arrangements  have  a  significant  impact  on  the  pressure  distribu¬ 
tion.  In  Case  5,  a  significant  amount  of  pressure  has  been  transferred 
to  the  middle  of  the  active  area  due  to  EQL  deformation. 

4.2.  About  numerical  solutions 

The  mesh  was  automatically  generated  by  Comsol  Multiphysics. 

Mesh  element  size  setting  ‘Normal’  was  used  for  Case  1  and  ‘Fine’ 
for  other  cases.  The  models  had  from  74  to  175  k  tetrahedral 
mesh  elements.  Conjugate  gradients  iterative  solver  with  geometric 
multigrid  preconditioner  was  used  and  solution  times  varied  from 
15  to  50  min  depending  on  the  model.  Relative  error  estimates  for 
nonlinear  solver  results  were  on  the  order  of  10-15  to  10-16.  To 
investigate  the  effect  of  mesh  resolution,  Case  5  was  solved  with  an 
alternative  mesh  consisting  of  323  k  elements.  However,  the  results 
were  identical  to  the  solution  with  the  coarser  mesh. 

4.3.  Experimental  results 

Measured  pressure  distributions  for  Cases  1-5  are  shown  in 
Fig.  5.  In  all  cases,  the  measured  pressure  distributions  show  the 
same  trends  as  the  model  predictions.  However,  there  are  notable 
quantitative  differences  in  areas  where  the  measured  pressure  is 
close  to  the  upper  or  lower  limit  of  the  pressure  sensitive  film’s 
range  (4-25  bar).  Below  4  bar,  the  pressure  does  not  register  on  the 
film,  and  the  colouration  reaches  saturation  at  25  bar.  Furthermore, 
data  averaging,  which  was  performed  to  reduce  noise,  obscures  the 
areas  of  highest  compression.  In  general,  all  measured  pressure 
distributions  indicate  lower  pressures  than  their  modelled  coun¬ 
terparts. 

4.4.  Model  validation  results 

The  experimental  results  in  Fig.  5  show  the  same  trends  as  model 
predictions  presented  in  Fig.  4.  For  numerical  analysis,  the  total 
forces  on  active  areas  were  determined  (see  Table  3).  Consistently, 
the  model  predicts  higher  compression  than  the  measurements 
indicate.  The  most  probable  cause  for  this  behaviour  is  the  inaccu¬ 
racy  of  material  parameters  used  in  the  model,  especially  Young’s 
modulus  for  the  current  collector  and  equalization  layer. 

The  correctness  of  the  shape  of  the  modelled  pressure  dis¬ 
tribution  was  evaluated  using  the  following  procedure.  First,  the 
pressure  for  both  the  modelled  and  measured  results  on  each  grid 
square  was  scaled  by  dividing  it  by  the  maximum  value  in  that  grid. 

Then  the  error  between  the  measured  and  modelled  was  calculated 
using  Eq.  (4).  Error  values  are  presented  in  Table  4.  - 

g> 

Error  =  100%  x  (measured  -  modelled)/modelled  (4)  .. 

The  largest  errors  in  each  case  are  seen  where  the  measured 
pressure  was  close  to  the  lower  end  of  the  measurable  range,  but 
where  the  measured  pressure  was  between  7  and  20  bar,  the  error 
is  generally  below  30%.  This  confirms  that  the  model  predictions 
are  qualitatively  correct.  Among  the  cases,  the  largest  errors  occur 

Table  3 

Total  compressive  force  on  active  area  for  measured  and  modelled  cases  (1<N). 


Table  4 

The  difference  between  the  measured  and  modelled  pressure  distributions.  The  grid 
coordinates  refer  to  those  in  Fig.  3. 


Grid 

coordinate 

1 

2 

3 

4 

5 

6 

4 

Case  1 
-51.5% 

-33.4% 

-16.0% 

-1.2% 

1.4% 

0.0% 

3 

-58.7% 

-54.5% 

-38.5% 

-23.2% 

-19.4% 

-5.7% 

2 

-60.1% 

-62.3% 

-62.6% 

-33.9% 

-20.9% 

-12.7% 

1 

-58.1% 

-58.0% 

-58.8% 

-33.7% 

-13.8% 

-16.2% 

4 

Case  2 

31.8% 

28.5% 

28.3% 

31.0% 

41.5% 

48.7% 

3 

40.7% 

25.0% 

19.4% 

11.2% 

24.6% 

57.4% 

2 

29.7% 

14.6% 

7.7% 

4.2% 

-0.4% 

43.4% 

1 

27.9% 

13.1% 

7.2% 

-7.0% 

-5.8% 

31.2% 

4 

Case  3 
7.3% 

-11.3% 

-6.8% 

0.7% 

23.4% 

14.5% 

3 

-3.2% 

-9.6% 

-9.4% 

3.8% 

20.9% 

5.0% 

2 

3.6% 

6.8% 

4.9% 

8.7% 

31.5% 

-3.3% 

1 

0.7% 

3.5% 

1.2% 

0.0% 

25.2% 

-12.1% 

4 

Case  4 

-25.7% 

-44.8% 

-37.6% 

-11.7% 

2.7% 

1.5% 

3 

-6.1% 

-8.2% 

-5.1% 

16.5% 

35.8% 

1.5% 

2 

11.5% 

14.0% 

9.3% 

15.4% 

42.8% 

-5.6% 

1 

5.6% 

5.9% 

-3.7% 

0.0% 

28.9% 

-7.2% 

4 

Case  5 
-21.6% 

-42.8% 

-42.0% 

-21.4% 

-8.8% 

-10.8% 

3 

-12.9% 

-12.0% 

-8.2% 

10.3% 

28.9% 

-9.2% 

2 

2.9% 

6.8% 

4.8% 

25.4% 

41.1% 

-16.9% 

1 

-11.4% 

-10.4% 

-7.8% 

15.4% 

35.5% 

-17.3% 

in  Case  1,  where  close  to  80%  of  the  compressive  force  was  exerted 
outside  the  active  area. 

4.5.  Sensitivity  analysis 

The  model’s  sensitivity  to  material  parameters  was  investigated 
by  solving  the  model  for  Case  5  using  various  parameter  values  and 
comparing  the  resulting  clamping  pressure  distribution  to  the  mea¬ 
sured  one.  Young’s  modulus  of  the  end  plate  (FEP),  equalization  layer 
(Feql)  and  current  collector  (FCc)  were  varied  one  at  a  time  from  25% 
to  175%  of  the  values  listed  in  Table  1.  The  resulting  distributions 


Fig.  6.  The  model’s  sensitivity  to  Young’s  moduli  of  the  end  plate,  equalization  layer 
and  current  collector.  The  moduli  are  given  in  percentage  of  values  in  Table  1. 
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were  analyzed  by  summing  the  squares  of  errors  (see  Eq.  (4))  for 
each  parameter  value.  The  parameter  value  with  the  corresponding 
minimum  value  of  the  squares  sum  produced  the  pressure  distri¬ 
bution  which  resembled  the  measured  one  the  most.  The  square 
sums  of  errors  are  presented  in  Fig.  6.  For  FEP  and  FEQE  the  mini¬ 
mum  squares  sum  was  found  at  100%  of  the  Table  1  values,  but  for 
Fee  the  minimum  was  at  50%  value.  This  indicates  that  the  Young’s 
modulus  used  for  the  current  collector  may  be  inaccurate.  Since  the 
error  was  relatively  small  and  a  more  reliable  value  was  not  avail¬ 
able,  the  estimate  for  Young’s  modulus  of  the  current  collector  was 
not  changed. 

5.  Conclusions 

A  FEM-model  which  predicts  the  compression  pressure  distribu¬ 
tion  and  component  deformation  in  the  outermost  cells  in  a  fuel  cell 
stack  has  been  developed  and  experimentally  validated.  It  shows 
good  qualitative  agreement  with  measured  pressure  distributions 
for  different  equalization  layers,  clamping  schemes  and  forces. 
Quantitatively,  the  model  predicted  higher  compression  pressures 
than  were  measured.  For  a  large  part,  this  is  due  to  the  limited  range 
of  the  pressure  sensitive  film  that  was  used.  Here,  it  was  observed 
that  the  colour  intensity  of  5-10%  of  the  pixels  in  the  scans  of  the 
pressure  sensitive  films  was  at  saturation  limit.  In  retrospect,  it  is 
clear  that  several  pressure  sensitive  films  with  different  pressure 
ranges  should  be  used  in  experiments  to  record  a  larger  portion  of 
the  pressure  spectra. 

Other  possible  sources  of  error  are  the  uncertainty  in  some  of  the 
material  parameters  used  in  the  model,  and  manufacturing  toler¬ 
ances  in  the  experimental  hardware.  It  has  been  observed  that  even 
small  errors  in  component  dimensions  can  cause  notable  variations 
in  the  pressure  distribution  [14,18].  A  preliminary  sensitivity  anal¬ 
ysis  was  performed  to  investigate  the  sensitivity  of  the  model  to 
material  parameters. 

In  the  future,  the  validation  process  will  be  developed  further, 
and  thermal  expansion  will  be  included  in  the  model.  This,  com¬ 
bined  with  accurate  material  property  data  library,  will  facilitate 
predicting  compression  pressure  distribution  reliably  even  at  ele¬ 
vated  temperatures. 
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